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1 Abstract 

The objective of this paper is to extend our recently developed highly parallelizable non- 
linear stable high order schemes for complex multiscale hydrodynamic applications to the 
viscous MHD equations. These schemes employed multiresolution wavelets as adaptive 
numerical dissipation controls to limit the amount of and to aid the selection and/or 
blending of the appropriate types of dissipation to be used. The new scheme is formulated 
for both the conservative and non-conservative form of the MHD equations in curvilinear 
grids. The four advantages of the present approach over existing MHD schemes reported 
in the open literature are as follows. First, the scheme is constructed for long-time inte- 
grations of shock/turbulence/combustion MHD flows. Available schemes are too diffusive 
for long-time integrations and/or turbulence/combustion problems. Second, unlike exist- 
ing schemes for the conservative MHD equations which suffer from ill-conditioned eigen- 
decompositions, the present scheme makes use of a well-conditioned eigen-decomposition 
obtained from a minor modification of the eigenvectors of the non-conservative MHD 
equations to solve the conservative form of the MHD equations. Third, this approach 
of using the non-conservative eigensystem when solving the conservative equations also 
works well in the context of standard shock-capturing schemes for the MHD equations. 
Fourth, a new approach to minimize the numerical error of the divergence-free magnetic 
condition for high order schemes is introduced. Numerical experiments with typical MHD 
model problems revealed the applicability of the newly developed schemes for the MHD 
equations. 


2 Motivation & Relevance 

Accurate numerical simulations of complex multiscale compressible viscous flows, espe- 
cially high speed turbulence combustion and acoustics, demand high order schemes with 
adaptive numerical dissipation controls. Standard high resolution shock- capturing meth- 
ods are too dissipative to capture the small scales and/or long-time wave propagations 
without extreme grid refinements and small time steps. An integrated approach for the 
control of numerical dissipation in high order schemes for the compressible Euler and 
Navier-Stokes equations has been developed and verified by the authors and collaborators 
[18, 19, 12, 20, 15]. These schemes are suitable for the problems in question. Basically, 
the scheme consists of sixth-order or higher non-dissipative spatial difference operators 
as the base scheme. To control the amount and types of numerical dissipation, an artifi- 
cial compression method (ACM) sensor or multiresolution wavelets are used as sensors to 
adaptively limit the amount and to aid in the selection and/or blending of the appropriate 
types of numerical dissipation to be used. This adaptive control of numerical dissipation 


2 


Bjorn Sjogreen and H. C. Yee 


is accomplished by a filter step after the completion of each full time step integration of 
the base scheme. 

MHD and/or plasma injections into a gas dynamics field could play an important role 
in drag reduction in highly maneuverable high speed combat aircraft. MHD plays a key role 
in space weather forecasting, and in the understanding of the dynamics of the evolution of 
our solar system and the main sequence stars, especially certain types of star formation. 
Although there exist a few well-studied second- and third-order high-resolution shock- 
capturing schemes for MHD in the open literature, these schemes are too diffusive and not 
practical for turbulence/combustion MHD flows. On the other hand, extension of higher 
than third-order high-resolution schemes to the MHD system of equations is not straight- 
forward. Unlike the hydrodynamic equations, the inviscid MHD system is non-strictly 
hyperbolic with non-convex fluxes. The wave structures and shock types are different 
from their hydrodynamic counterparts. Many of the non-traditional hydrodynamic shocks 
are not fully understood. Consequently, reliable and highly accurate numerical schemes for 
multiscale MHD equations pose a great challenge to algorithm development. In addition, 
controlling the numerical error of the divergence- free condition of the magnetic field for 
high order methods has been a stumbling block. Lower order methods are not practical for 
the astrophysical problems in question. We propose to extend our hydrodynamics schemes 
to the MHD equations, gaining several desirable advantages over commonly used MHD 
schemes. 

The present paper is Part I of a series of papers on the subject. Here, we present our 
new scheme with typical ideal MHD test cases to validate our approach and to illustrate 
the representative performance of the scheme. Since there are many variants in the min- 
imization of the divergence of the magnetic field (V- B) numerical error, Part II of the 
companion paper [21] is devoted to extensive comparison of these variants. In yet another 
forthcoming paper, Part III, complex multiscale high speed turbulent astrophysical appli- 
cations will be sought. Throughout the paper, the term “V*B numerical error” refers to 
the “amount of non-zero value of the discretized form of V • B of the underlying scheme” . 


3 New Approach 

Aside from possessing the desirable property of adaptive numerical dissipation control of 
our proposed MHD scheme, the new approach consists of several additional new concepts. 
First, our filter scheme employs high order central spatial difference operators that preserve 
the divergence-free condition of the magnetic field as the base scheme. The only possible 
violation of the divergence-free condition is from the filter step. Second, a new idea in 
applying our filter procedure for minimizing the numerical error of V- B for high order 
schemes is proposed. For certain flows, we are able to preserve the divergence-free condition 
for the complete scheme. The third new concept is as follows: 

The MHD equations are a system of non-strictly hyperbolic conservation laws. The 
non-convexity of the inviscid flux vector results in corresponding Jacobian matrices with 
undesirable properties. It has previously been shown by Powell et al. [10] that an “al- 
most” equivalent MHD system in non-conservative form can be derived. In order to have 
a better conditioned eigensystem for the application of high-resolution shock-capturing 
schemes, they adjoined ad hoc non-conservative terms to the conservative equation. Ac- 
tually, the MHD equations can be derived from basic principles in either conservative or 
non-conservative form (Marcel Vinokur of NASA Ames, private communication, 1996). 
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We formulate our new scheme together with the Cargo &: Gallice [3] form of the MHD 
Roe approximate Riemann solver in curvilinear grids for both versions of the MHD equa- 
tions. A novel feature of our new method is that the well-conditioned eigen-decomposition 
of the non-conservative MHD equations, with a minor modification, is used to solve the 
conservative equations. This new feature of the method provides well-conditioned eigen- 
vectors for the conservative formulation, so that correct wave speeds for discontinuities are 
assured. Another advantage of using the non-conservative eigen-decomposition to solve the 
conservative equations is that our scheme exhibits smaller V- B numerical error. It will be 
shown that this approach, using the non- conservative eigensystem when solving the con- 
servative equations, also works well in the context of standard shock-capturing schemes 
involving the use of the eigen-structure of the MHD equations. 


4 Description of the Scheme 

A full description of our adaptive low dissipative high order scheme for the Euler and 
Navier-Stokes equations can be found in [18, 19, 12, 20, 15]. Here, we describe the extension 
of this scheme to the MHD equations with the blending of high order nonlinear filters 
and high order linear filters. An important ingredient in our method is the use of the 
dissipative portion of high-resolution shock-capturing schemes as nonlinear filters. These 
nonlinear filters involve the use of approximate Riemann solvers. We will therefore first 
present a new form of high-resolution shock-capturing schemes for the conservative MHD 
equations using the non-conservative eigensystem. 


4.1 Preliminary 

Consider the 3-D non-conservative form of the ideal MHD equations 
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where the velocity vector is denoted u = (u,v i w) T , the magnetic field vector is B = 
(. B x , By , B z ) t , p is the density, and e is the total energy. The vector on the right hand side 
is the non-conservative portion of the MHD equations. The notation B 2 = B 2 + By + B 2 
is used. The pressure is related to the other variables by 

P = (7 - l)(e - \p{u 2 + v 2 + w 2 ) Bl + B 2 )). 

For plasmas, 7 is usually equal to 1.667 (for Mon-atomic gases). 

We will also consider the conservative MHD equations, obtained by setting the right 
hand side equal to zero. The non-conservative term is proportional to the divergence of 
the magnetic field (V* B). Physically, it is zero if V- B = 0 initially. 

In symbolic form, the non-conservative form can be written as 


4 


Bjorn Sjogreen and H. C. Yee 


/7 £ + V*E = S, 

where U is the corresponding state vector, F is the conservative inviscid flux vector and 
S is the non-conservative portion of the equations in (1). 

There has been much concern in preserving the divergence- free condition throughout 
the computation, see e.g., [5, 16] and references cited therein. The V*B — 0 condition is 
an initial constraint for the MHD equations. The divergence- free condition is not part of 
the MHD differential system, and it is very different from the divergence-free condition 
for the incompressible Euler or Navier-Stokes equations. The divergence-free condition for 
the incompressible Euler or Navier-Stokes is part of the differential system. It is needed to 
close the system and must be enforced explicitly. For the MHD equations, as long as the 
magnitude of the discretized form of V* B is on the order of truncation error, and goes to 
zero when the grid is refined, there should be no problem. Unfortunately, not all schemes 
preserve the divergence- free condition. This is especially true for standard upwind and 
high-resolution shock-capturing schemes that involve some form of Riemann solver. 

When using pure centered difference operators, it is trivial to see that the divergence of 
B is perfectly preserved. Take for example the semi-discrete approximation of the magnetic 
field equations 

dB x (t)^ j ■ 
dt 

B X ij t k v iJ,k) + Dk(B zi j k Uij'k “ BxiJ’kWiJ^k) ” 0 
dt 

j — B Vi - k Uij t k) + Bk(B zi j k Vi 7 j^k — B Vi . k Wij 7 k) = 0 

dB g (t) iti , k 
dt "*■ 

Di(B X ij ik Wij,k B Z i i j i k u ij, k ) 4- Dj(By iJ k Wi J)k — B Z ij k Vij 7k ) = 0, 

where Di , Dj and D k denote finite difference operators acting in the i-, j-, ^-directions, 
respectively. 


Forming the divergence by taking the sum of Di on the B x equation, Dj on the B y 
equation, and D k on the B z equation gives 


d(DiB x + DjB y 4- D k B z )i t j y k __ ^ n D _ inno 

— DiDjBy^m^ 4 - DjDiB yi j k Uij )k + ... — 0 , 


where the dots denote several additional terms of the similar form as the first two. All 
these terms disappear, since the difference operators along different coordinate directions 
commute, i.e., DiDj = DjDi. 

When the solution is smooth, we can use a high order centered difference operator and 
perfectly preserve the divergence-free condition. In this case the result will be the same, 
whether we solve the conservative or non-conservative equations. 


From now on, the discussion pertains to schemes involving the use of Riemann solvers or 
the eigen-structure of the MHD equations. For convenience of presentation, we will describe 
our numerical methods for 1-D problems on a uniform grid. The 1-D MHD equations 
become 

U t + F(U) X = 0 (conservative) 

U t + F(U) X = S (non-conservative) 
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Let A(U) denote the Jacobian dF/dU with the understanding that the present U, F and 
S are the 1-D counterpart of the 3-D description above. For later discussion, we write the 
non-conservative terms S = N(U)U X , 

Presently, there are basically two camps in solving the MHD equations; namely, solving 
the conservative or more recently, popularized by [10], the non-conservative form. In both 
cases high-resolution shock capturing methods suffer from the need to perform extra work 
to keep the V* B numerical error to machine zero. Three of the popular approaches are to 
minimize the V-B numerical error by augmenting an extra PDE to the system [5], using 
the staggered approach of K.S. Yee [22, 6, 7], and a projection method [11]. There is, 
however, a key advantage to solving the conservative equations over the non-conservative 
equations since the conservative form guarantees correct propagation speeds and locations 
of discontinuities. 

It can be shown that the “non-conservative Jacobian” A(U) — N(U) has real eigen- 
values, and a complete set of eigenvectors [3]. The conservative Jacobian A(U) has real 
eigenvalues, but there exist states U for which A(U) does not have a complete set of 
eigenvectors. Therefore, defining an upwind scheme based on A(U) is problematic. 

In the next subsection, we show a simple way of constructing upwind or high-resolution 
shock-capturing schemes that use approximate Riemann solvers for the conservative equa- 
tions without the use of its ill-conditioned conservative eigensystem. This will set the stage 
for a conservative filter approach involving the use of the resulting nonlinear dissipative 
portion of these high-resolution shock-capturing schemes. 

4.2 Conservative and Non- Conservative Formulations Involving the Use of 
Approximate Riemann Solvers 

For simplicity of presentation, we illustrate the idea for a first-order upwind scheme in- 
volving the use of the MHD counterpart of Roe’s approximate Riemann solver for the 
gas dynamic equations. Here we show that the construction of an upwind scheme for the 
conservative equations is possible by using the non-conservative eigen- decomposition. Fur- 
thermore, it will be shown by numerical experiments that the conservative method will 
give smaller generation of spurious divergence of the magnetic field, than that obtained 
when solving the non-conservative equations. 

In Cargo & Gallice [3] a matrix, A j+ i / 2 = A(Uj, Uj+i), consistent with the Jacobian, 
and having the “property U” of Roe for the MHD equations 

^i+i - Fj = Aj+i^iUj+i - Uj) 

is constructed. Unlike the case of standard gas dynamics, this matrix cannot be written 
as the Jacobian at an average state of A. By inspecting terms in the matrix, they use 
the decomposition A J+1 / 2 = (A, +1 / 2 — Nj+ 1 / 2 ) + 1 / 2 , where iVj+ 1/2 is consistent 

with N(U). This is due to the fact that there is a simple relationship between the eigen- 
structure of the conservative and the non-conservative forms. Here N(U) is part of the 
non-conservative terms S = N(U)U X . 

It turns out that seven of the eigenvalues and eigenvectors are identical for the Jacobian 
matrices {A J+1 / 2 ) and (Aj+i/ 2 - N j+l / 2 ) [8]. Assuming the eigenvalues are arranged in the 
increasing order, the difference between the two matrices lies only in the eighth eigenvector. 
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Furthermore, the eighth eigenvector of the conservative system can sometimes coincide 
with one of the eigenvectors, thereby, making it impossible to define the MHD Roe’s 
approximate Riemann solver in the standard way. The eigenvectors of the non-conservative 
Jacobian (A j+x / 2 — Nj +x / 2 ) are very simple. 

We proposed to use the non-conservative form but with the fourth eigenvalue (the 
degenerate eigenvalue) replaced by zero for the conservative form, since only the eighth 
eigenvector of the non- conservative form is not the same as the eighth eigenvector for the 
conservative form. The incorrect eigenvector will be multiplied by an eigenvalue which is 
close to zero (the eigenvalue will not be exactly zero when an entropy correction is used; 
see below). Thus the effect of a “false” eigenvector will be small. By using the eighth 
eigenvector of the non-conservative system instead, the difficulty will be avoided. 

Let k — 1 ,...,8 denote the eigenvalues of the conservative Jacobian Aj+ x / 2 . 

Let A j + i / 2 an( i ijL}_ i/2 denote the eigenvalues and eigenvectors, respectively, of the non- 
conservative Jacobian Aj +X f 2 — Nj+ x / 2 . With the above justification, we define the nu- 
merical upwind flux of the conservative equations as 

1 1 A 

*7+1/2 — +Fj) — x 

z 1 i 

where a^ +1 / 2 is the jump of the respective characteristic variables of Uj+ X —Uj in the eigen- 
vector basis, and q(x , e) is the entropy-corrected absolute value, e.g., q(x , e) = max(|a;[, e), 
where e is a small parameter. Here fj +1 / 2 i s the minor modification of non-conservative 
eigenvectors discussed previously (only the eighth eigenvector is modified). 

For the non-conservative equations, we define the corresponding numerical flux 

1/2 = \(F j+ 1 +Fj) - i <l( X j+l/2’ e ) a j+l/2 T j+l/2- 

z k= 1 


Making use of the identity Aj + X / 2 = (Aj + i / 2 — A^+1/2) + Nj+i/ 2 , the semi-discrete non- 
conservative upwind approximation of the MHD equations can be written as 

^ = -^(^ 1/2 - 1 / 2 ) + ^(N j+ 1 / 2 (U j+1 - U S ) + N^Uj - U^)). ( 2 ) 

The semi-discrete conservative upwind approximation of the MHD equations is the usual 
form 

f-sU'w-W < 3 > 

The conservative and non-conservative TVD and WENO schemes used in this paper are 
generalized from either of these two first-order schemes in a standard way. 

4.3 High Order ACM and Wavelet NonLinear/Linear Filter Method 

Our high order filter method consists of two stages, a divergence-free preserving base 
scheme stage and a filter stage which can be divergence-free preserving depending on the 
type of filter operator being used and the method of applying the filter step. 



High Order Methods for MHD 


7 


Divergence- Free Preserving Base Schemes: The first stage of the numerical method 
consists of a time step by a spatially high order non-dissipative and high order temporally 
base scheme operator L (e.g., a divergence-free preserving sixth-order central in space and 
high order linear-multistep or fourth-order Runge-Kutta in time), 

U* = L(U n ), 

where U n is the numerical solution vector at time level n. When necessary, a high order lin- 
ear numerical dissipation operator can be used. For example, a divergence-free preserving 
eighth-order linear dissipation with the sixth-order centered base scheme to approximate 
F(U)x is written, in the 1-D case, as 

8F 

— « -DoeFj - dAx 7 (D + D-) 4 U j , (4) 

where Dqq is the standard sixth-order accurate centered difference operator, and D+D- is 
the standard second-order accurate centered approximation of the second derivative, d is 
a tunable parameter. The operators are modified at boundaries in a stable way [12]. (Note 
that for more than two level linear-multistep methods, the L operator involves the corre- 
sponding number of time levels.) This highly accurate base scheme is constructed to numer- 
ically preserve the divergence-free condition of the magnetic field (to the level of round-off 
error) for curvilinear grids. The only possible source of violation of the divergence-free con- 
dition is from the filter step (assuming the correct handling of the physical and numerical 
boundary conditions). 

Adaptive Numerical Dissipation Filters: After the completion of a full time step of 
the divergence-free preserving base scheme, the second stage is to filter the solution by 
the dissipative portion of a high-resolution shock-capturing scheme. The filter operator 
can be written as (assume 1-D for ease of illustration) 

The filter numerical flux vector is 

Hj+i/2 = ^j + l/2^ + 1/2* 

Here A/+ 1/2 is the matrix of right eigenvectors of the Jacobian of the non-conservative 
MHD flux vector (A, +1 / 2 — iVj+ 1 / 2 ) as discussed in the previous subsection. The Hj + i/ 2 
are also evaluated from the same characteristic quantities derived from these eigenvectors. 
Denote the elements of the vector Zfj +1 / 2 by h j+1 f 2 , l = 1,2, ..., 8. They have the form 

hj+1/2 — (^^)j+l/2 (^'+1/2) “ ( sL )j+i/2^+i/2* 

Here (s N ) l j+i/ 2 and (s L ) l j+ 1 / 2 are sensors to activate the higher order nonlinear filter and 
linear filter, respectively. For example, (s JV )J +1 / 2 i s designed to be zero in regions of smooth 
flow and near one in region with discontinuities. (s N )* +1 / 2 varies from one grid point to 
another and is obtained either from a wavelet analysis of the solution (WAV-filter scheme), 
or from a gradient-based detector (ACM-filter scheme) [18, 19, 12, 20, 15]. We have in the 
previous numerical experiments [20] used (s L ) l J+1 / 2 = 1 ~ (s N ) l j+ 1/2 > but other choices are 
possible. See [20] for a discussion. The functions <^ +1 / 2 an d dj+ 1/2 are the nonlinear and 
linear filters for the Zth-characteristic wave. 
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The nonlinear filter <j> l j+1 , 2 = <^ +1 / 2 — ^+ 1/2 is dissipative portion of a uniformly 
high order high-resolution shock-capturing scheme for the /th-characteristic wave. Here 
g l j + 1/2 and ^ +1 / 2 are numerical fluxes of the uniformly high order high-resolution scheme 
and a high order central scheme for the Zth characteristic, respectively. It is noted that 
b l j+l/2 might not be unique since there is more than one way of obtaining F° r the 

forms of the <^ +1 / 2 nsed in the numerical experiment section, see [18, 19, 12, 20, 15]. For 
example, the form of Harten and Yee and symmetric TVD schemes are already in the 
proper form in the sense that they are written in a central differencing portion ^ +1 / 2 and 
a nonlinear dissipation portion 0j +1 / 2 - No work is required to obtain <£j + 1 / 2 in this case. 

The linear filter d*. +1 ^ 2 consists of a tuning parameter df. For a sixth-order spatial 
base scheme, the eighth-order central dissipation for the linear filter is used. For the test 
cases to be shown later, the majority of the computations only use the nonlinear filter. It 
is not necessary to use the linear filter, i.e., setting (s L ) l j+l / 2 = 0- 

When Runge-Kutta time stepping is used, the filter is usually applied after the com- 
pletion of each full time step, or for extremely high speed flows, applied after each stage 
in the Runge-Kutta method. In all of the numerical examples shown later, we filter the 
solutions after the completion of each full Runge-Kutta time step. 

Options in Filtering the Magnetic Field Equations: In order to minimize the nu- 
merical error of the divergence-free magnetic condition, the nonlinear filter step only acts 
on the hydrodynamic portion of the equations. To further improve nonlinear stability 
and accuracy, if necessary, we might employ high order linear filters that preserve the 
divergence-free property for the magnetic field portion of the equations. For problems 
containing moderate to strong shocks that are magnetically dominated, it might be nec- 
essary to apply a nonlinear filter on the magnetic field equations as well. In this case, the 
scheme, through wavelet sensors, still maintains the divergence- free preserving property in 
smooth regions of the flow. Construction of accurate schemes for magnetically dominated 
shocks is a subject of ongoing research. The next section illustrates the representative 
property of our scheme for mainly shock wave-dominated MHD test cases. 


5 Shock Wave Dominated MHD Numerical Examples 

For illustration purposes, numerical examples using sixth-order central spatial and fourth- 
order temporal discretizations as the base scheme with adaptive numerical dissipation con- 
trols for four typical MHD model problems are included. These are shock wave-dominated 
examples. The purpose is to validate our new scheme for conventional MHD model prob- 
lems. Since most of these examples are not turbulence/ combustion flows, the full capability 
of the new scheme is not utilized. Thus, we do not expect the new scheme to exhibit no- 
ticeable improvement in stability and shock-resolution over conventional schemes. In a 
forthcoming paper, multiscale MHD examples will be included and new features of the 
present scheme will surface. 

The sixth-order base scheme together with the nonlinear /linear filter with wavelet 
sensor will be denoted WAV66. When a more conventional gradient based sensor is used, 
the scheme is denoted ACM66. If high order linear numerical dissipation is also used in the 
base scheme, the methods will be denoted WAV66+AD8 and ACM66-f AD8 respectively. 
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The strength of the eighth-order dissipation will be denoted by a tunable coefficient d , as 
in (4). When a high order filter is used, its dissipation strength will be denoted b yd/. 

For comparison, the fifth-order weighted ENO scheme (denoted WEN05), described in 
[9] will also be used. Classical fourth-order Runge-Kutta time stepping is used for all sixth- 
order schemes, as well as for the WEN05 scheme. The second-order TVD schemes Harten- 
Yee and MUSCL are integrated in time by the second-order TVD Runge-Kutta method. 
The conservative (3) and non-conservative (2) versions of the TVD schemes have both 
been implemented and will be compared below. All the methods use the same approximate 
Riemann solver of Cargo & Gallice. Unless indicated, all of the computations use uniform 
grids and solve the conservative form of the MHD equations with the non-conservative 
eigenvector decomposition as discussed in Section 4.2, and Section 4.3 for the filter scheme 
(ACM66 and WAV66). In addition, for ACM66, ACM66+AD8, WAV66 and WAV66+AD8, 
unless indicated, the nonlinear filter are applied on the full MHD system on numerical 
solutions shown. 

The V* B numerical error is obtained by approximating the spatial derivatives by sixth- 
order centered differences for WAV66, ACM66 and WEN05, whereas the corresponding 
V- B numerical error is obtained by second-order centered differences for the second-order 
TVD schemes (MUSCL and Harten-Yee). 


5.1 A 1-D Riemann Problem 

Here we consider the same 1-D Riemann (shock-tube) problem as examined by Brio & Wu 
[2]. This Riemann problem is chosen such that the x and 2 components of the magnetic 
field, B x and B z , are constant throughout the time evolution. In other words, the magnetic 
field does not change direction, only magnitude. The initial data for the left (L) and right 
(R) states are 
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where p is the density, u, v and w are the velocity components, and e is the total energy per 
unit volume. The B XJ B y and B z are the three magnetic field components. The density and 
the y-component of the magnetic field of a reference solution obtained by a second-order 
Harten-Yee TVD scheme for the non-conservative MHD form of the equations, using 16000 
uniform grid points, are shown in Fig. 1. Note that the eighth wave is not present and the 
rotational discontinuities are also not present in this 1-D co-planar Riemann problem. We 
are thus left with 5 waves. 

Comparisons of MUSCL, Harten-Yee and WEN05 with the present scheme were con- 
ducted to validate our approach. All computations were run on the domain 0 < x < 1 to 
time 0.12, using a fixed CFL number of 0.2 and using the same non-conservative form of 
the MHD equations. A uniform grid of 400 points is used for the comparison. The value 
of the gas law parameter 7 was 1.6667. Figure 2 shows the density using the MUSCL 
second-order TVD scheme with the van Albada limiter (solid: reference solution; circle: 
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Fig. 1. Reference Solution for the 1-D Riemann problem. Density and y-component of the mag- 
netic field computed by Hart en- Yee at time 0.12 using a uniform 16000 grid. 
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Fig. 2. 1-D Riemann problem. Density computed by MUSCL with van Albada limiter using a 
uniform 400 grid. 
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Fig. 3. 1-D Riemann problem. Density computed by WEN05 using a uniform 400 grid. 


MUSCL). The upper left plot shows the solution on the entire computational domain. The 
sub-figures zoom in on different key parts of the domain. 

Figure 3 shows the computation by WEN05. Figure 4 shows the computation using 
the WAV66+AD8 filter scheme (sixth-order spatial base scheme, second-order nonlin- 
ear filter and blending with a small amount of high order linear filter). It can be seen 
that the WAV66-f AD8 exhibits a slightly higher accuracy than the MUSCL and WEN05 
schemes. Note that for this rapidly developing short time integration flow without any 
structure away from discontinuities, it is expected that any good second-order shock- 
capturing scheme would perform well. The small over- and undershoots are characteristic 
of all studied schemes (and numerical results published in the open literature) even with 
the use of 8000 uniform grid points. Adaptive grid refinement or the use of extremely fine 
grids (e.g., 16000 grid points) is necessary to minimize the small oscillations. 

5.2 Orszag-Tang Vortex 

The purpose of this example is to illustrate the advantage of the filter approach over the 
standard shock-capturing scheme in minimizing the numerical error of the divergence-free 
magnetic condition. The advantage of the present conservative scheme over the conven- 
tional conservative scheme is also illustrated. 

The 2-D Orszag-Tang vortex problem consists of periodic boundary conditions with 
smooth initial data 
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WAV66, van Albada e=0.2, ds=Q.0001, df=0.002 
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Fig. 4. 1-D Riemann problem. Density computed by WAV66+AD8 with van Albada limiter using 
a uniform 400 grid. 
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The computational domain is 0 < x < 2i r, 0 < y < 2i r. The computation stops at time 
T= 3.14 7 r) using a CFL number of 0.6. 

Density contours of the solution for three levels of grid refinement (201 x 201, 401 x 401, 
and 801 x 801) by the Harten-Yee scheme solving the non-conservative MHD equations 
are displayed in Fig. 5. A grid convergence study of the WEN05 scheme, integrated in 
time with the classical fourth-order accurate Runge-Kutta method is shown in Fig. 6. The 
Harten-Yee scheme and the WEN05 schemes do not perfectly conserve the divergence of 
B. Here B is the 3-D magnetic field vector. We expect a higher order method to give 
a smaller error in the V-B. The development of the L 2 norm of V-B in time for the 
Harten-Yee scheme and the WEN05 scheme using 801 x 801 uniform grid points is shown 
in Fig. 7. As long as the solution is smooth, it is clear that a more accurate method gives 
a smaller increase in V*B. However, when shocks have formed, the error is due to the 
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Harten-Yee, 201x 201 Harten-Yee, 401x 401 



Fig. 5. Density contours: Grid refinement of the Orszag-Tang problem by Harten-Yee 


amount of the nonlinear dissipation of the scheme in the vicinity of discontinuities. This 
is clearly seen in Fig. 7. 

The L 2 norm of the numerical V B as a function of time is shown in Fig. 9 for the 
three levels of grid refinement. It is clear that the norm of the divergence of the magnetic 
field does not converge to zero as the grid is refined. The errors in V B are associated 
with shock formation, and should not be expected to diminish, since, unlike our filter 
approach, the Harten-Yee scheme is applied to the entire MHD equations. Figure 8 shows 
the computation by WAV66+AD8 using uniform grids of sizes 101 x 101, 201 x 201, and 
401 x 401. 

Figure 10 shows the L 2 norm of V* B for the Harten-Yee scheme (dashed curve) and 
the WAV66 scheme (solid curve) for the non-conservative form of the equations. Perfect 
preservation of the divergence-free condition is obtained by the WAV66 before shock waves 
develop. This is clearly visible in Fig. 10. Up to time around 0.7, the solution is smooth 
enough that the nonlinear filter is not very active. The norm of the V* B numerical error 
is then on the level of round-off error. This is due to the fact that, although the spatially 
sixth-order central base scheme of WAV66 preserves the divergence-free condition, the 
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WAV66+AD8, 101x101 



WAV66+AD8, 201x 201 



x 


WAV66+AD8, 401x 401 



x 


Fig. 8. Grid refinement study for the Orszag-Tang problem. Density contours by WAV66-I-AD8, 
d = 0.005. 



Fig. 9. Orszag-Tang problem. L 2 norm of numerical V* B vs. time on grids of increasing refine- 
ment by Harten-Yee. 
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Fig. 10. Orszag-Tang problem. L 2 norm of numerical V* B vs. time, by WAV66 and Harten-Yee. 



Fig. 11. Orszag-Tang problem. L 2 norm of numerical V-B vs. time. Conservative and non- 
conservative MHD by Harten-Yee. 


nonlinear filter, even though not too active, introduces a numerical error in the V-B of 
the same size as the truncation error of the approximation. At later times the nonlinear 
filter dissipation of WAV66 quickly takes the V* B numerical value to that of the Harten- 
Yee scheme after the onset of shock waves. 

A comparison of numerical results using MUSCL, ACM66 and ACM66-I-AD8 were 
also conducted. Since the accuracy of the numerical solution by the MUSCL method was 
inferior to the rest of the methods (as for all test cases), we will not discuss the MUSCL 
scheme from here on. The ACM66 and ACM66+AD8 are slightly more stable but exhibit 
a slightly larger V-B numerical error than WAV66 (figures not shown). 

The development of the norm of the numerical V- B in time is shown in Fig. 11 for the 
conservative and the non-conservative forms of the MHD equations. The computations use 
the same Harten-Yee scheme and the same eigen-decomposition of the non-conservative 
equations to solve both forms of the MHD equations. The conservative form leads to a 
smaller V-B numerical error. The same result, that the conservative form of the MHD 
system gives smaller spurious V-B, has been observed for all test problems we have solved. 

From the comparison in Fig. 10, we conclude that it might be possible to preserve 
the divergence-free condition if we avoid applying the filter on the three components of 
the magnetic field equations. This requires that the filtering of the remaining components 
be enough to suppress spurious oscillations. The linear high order dissipation of the base 
scheme is of constant strength, and therefore does not affect the numerical divergence of 
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the magnetic field. A proper choice of the strength of this high order dissipation helps 
keep spurious oscillations low when switching off the nonlinear filter on the magnetic field 
equations. (Here after refers to “without nonlinear filter on B”.) In all of the numerical 
experiments (except the 1-D Riemann problem), linear filters are not activated on all of 
the field equations. We show results of this technique in Fig. 12, where the norm of the 
divergence of the magnetic field is plotted vs. time for the WAV 66+AD8 method with (solid 
curves) and without (dashed curves) applying the nonlinear filter on the three magnetic 
field components of the MHD system. The dissipation was d = 0.005 for the fully filtered 
computation, and d = 0.001 for the computation without nonlinear filter on the magnetic 
field. Results from the three different uniform grids 101 x 101, 201 x 201, and 401 x 401 are 
shown for each method. The solutions at 401 x 401 grid points are shown for the methods 
with and without nonlinear filter on the magnetic field in Fig. 13. 



Fig. 12. Orszag-Tang problem. Divergence B versus time for uniform grids of 101 x 101, 201 x 201, 
and 401 x 401 with nonlinear filter on the B field equations (solid), without nonlinear filter on 
the B field equations (dashed). 


The solution without nonlinear filter on B has fewer oscillations than the filtered 
solution. This seems contradictory, but is due partially to the fact that the base scheme for 
the non-filtered B computation uses a different strength of an eighth-order linear numerical 
dissipation. More importantly, not filtering B preserves the V-B property (Fig. 12). In 
our companion paper Part II [21], we will illustrate that no filter on B might not always 
work for all problems. Sometimes we need to filter B to maintain nonlinear stability. 

5.3 A 2-D Riemann Problem 

The 2-D Riemann problem consists of four constant states at time zero, as shown in Fig. 14. 
The values of the initial states were taken from [5], and are given by 

(p pu pv pw e B x By B z )++ = 

(0.9308 1.4557 -0.4633 0.0575 5.0838 0.3501 0.9830 0.3050) x > 0 y > 0 
{p pu pv pw e B x By B z )_+ = 

(1.0304 1.5774 -1.0455 -0.1016 5.7813 0.3501 0.5078 0.1576) x < 0 y>0 
(p pu pv pw e B x By B z ) — = 

(1.0000 1.75 - 1.0000 0.0000 6.0000 0.5462 0.5078 0.2539) x<0y<0 
(p pu pv pw e B x By B z )+ _ = 

(1.8887 0.2334 - 1.7422 0.0733 12.999 0.5462 0.9830 0.4915) x > 0 y < 0 
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WAV66+AD8, 401x 401 


WAV66+AD8, 401x 401 . No B-filter 



Fig. 13. Orszag-Tang problem. Density contours by WAV66+AD8 with filter on the B field 
equations (left) and solution without filter on the B field equations (right) on magnetic field 
components. 


The boundaries are treated as open boundaries. The problem is solved on the domain 


2D Rt*rrwv> probtem. Four states at tkns z*ro 


u 

u 



u_ 

■ 




-t* * 1 1 

-1 -0.5 0 0.5 1 


Fig. 14. Schematic of the initial data for the 2-D Riemann problem. 


— 1 < x < 1, — 1 < y < 1 to time T= 0.2. The purpose of this example is to show 
the advantage of solving the conservative equations on high-resolution shock-capturing 
methods for a more than one-dimensional Riemann problem. 

Density contours on grids of increasing refinement, computed by the second-order 
Harten-Yee TVD scheme applied to the conservative MHD equations are given in Fig. 15. 
The same computation solving the non-conservative MHD equations is shown in Fig. 16. 

From the contours plots, there appear to have very little difference between the two 
methods. However, the L 2 -norm of V-B vs. time shown in Fig. 17 illustrates that the 
solution of conservative equations has a smaller numerical V- B than the non-conservative 
equations. A contour plot of V-B is displayed in Fig. 18. The contours are placed at 
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Fig. 15. Density contours of the 2-D Riemann problem. Grid refinement by Harten-Yee solving 
the conservative MHD equations. 


Hart»n-Y»«, 201x201 , Noo-conaarvattv* Haiten-YM, 401x401 , rvabv* Hart*n-Y*«, SOIxflOl , Non-oooaarvaHva 



Fig. 16. Density contours of the 2-D Riemann problem. Grid refinement by Harten-Yee TVD 
solving the non-conservative MHD equations. 



Fig. 17. 2-D Riemann problem. L 2 norm of V B vs. time. Harten-Yee using a 201 x 201 grid 
solving both the conservative and non-conservative MHD equations. 
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the same levels in the two plots. As expected, both methods have V • B numerical errors 
concentrated at the discontinuities. The better performance of the conservative method in 
this respect seems to be due to better handling of generation of V*B at discontinuities. 




(a) Conservative (b) Non-conservative 

Fig. 18. 2-D Riemann problem. V*B contours by Harten-Yee solving both the conservative and 
non-conservative MHD for a 201 x 201 uniform grid. 


5.4 Vortex Pairing 

Vortex pairing in a mixing layer was simulated in [18], where the compressible Navier- 
Stokes equations were solved. The set-up of the problem is given in Fig. 19. The Navier- 


Vortex Pairing in a Time-Developing Mixing Layer 

(Me=0.8, Re=1000, T^OOK, Prandtl #=0.72) 


Normalized with vortlclty thickness: 

r = u i ~ u * 

(< du/dy )^ a 

T & c: determined by assuming constant stagnation enthalpy 
!£■ 

u 1= 1, u 2 =-1,T 1 =T 2i 

Initai shear profile: U = 0.5 tanh (2y) 

Crocco-Busemann: 1 

Initial perturbations: 

2 

v> ~ S °* co»(2tri*/X. + 


BC : Periodic in x, slip walls in y 
2^32* X y 3inh(6yT/) 

uniform in z 


Xy = 100, by = 3.4 



X. =30,5=10 
oa = 0.05, = — t/2 

= 0.01, fa = — jr/2 



Fig. 19. Schematic of the vortex pairing problem. 


Stokes equations are solved with a Reynolds number of 1000. The grid is mildly stretched 
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towards the line y = 0 with a sine-hyp stretching function. The boundary conditions axe 
periodic in the x-direction and open on upper and lower ^/-boundaries. 

As seen in Fig. 19, the initial flow consists of two constant states separated by a mixing 
layer. A small perturbation, which generates two vortices, is introduced at time zero. The 
two vortices are approaching each other. The problem is solved to a time when the vortices 
have merged. 

Here, generalization of the hydrodynamic problem to MHD by imposing an initially 
constant magnetic field in the x-direction, B x = 0.1, is used as a test case. In addition, 
we are solving the non-ideal MHD equations in the sense that a resistivity coefficient of 
a — 100 and the same Reynolds number (1000) and Prandtl number (0.72) that were used 
in the non-magnetic case [18] are used here. The non-ideal MHD terms are not shown 
in (1). See [1] for the form. The purpose of this numerical experiment is to show the 
performance of the scheme for an non-ideal MHD. 

The viscous terms are discretized by the sixth-order central differencing for WAV66, 
ACM66 and WEN05, and second-order central differencing for MUSCL and Harten and 
Yee. 

We solve the problem to time 90. Figure 20 shows temperature contours at time 90 of 
the solution, computed by the Harten- Yee second-order TVD scheme. The grid is refined 
from 201 x 201 grid points to 1601 x 1601 points. Grid convergence is reached without 
problems, due to a fairly small Reynolds number and conductivity. 

The advantage of a high order method to resolve the scales is seen in Fig. 21, where 
a grid convergence study for the sixth-order WAV664-AD8 method is shown. The dissi- 
pations d = 0.001 and the CFL number 0.6 were used. By comparing Figs. 20 and 21, 
we conclude that the sixth-order scheme captures all features of the flow correctly with 
201 x 201 grid points. The solutions (figures not shown) using ACM66 (with or with- 
out AD8) is similar to WAV66+AD8 except ACM66 is more stable than WAV66+AD8. 
WAV66 (without AD8) is not stable. The cause of the difference is very complex and is 
currently under investigation since the non-ideal MHD is an incompletely parabolic PDE 
system. The findings will be reported in a future publication. The second-order scheme 
needs to use 401 x 401 points to obtain a similar quality of the solution. This is in agree- 
ment with similar investigations for non- MHD computations [13]. Finally, in Fig. 22, we 
show a comparison of the development of L 2 norm of V-B in time for the second-order 
accurate Harten- Yee scheme, and the sixth-order WAV66+AD8 scheme. Computations 
using the WENOS schemes for this problem require an excessively small time step and are 
extremely slow to converge (solution not shown). 


6 Concluding Remarks 

The generalization of our low dissipative high order schemes for hydrodynamics equations 
to MHD equations has been presented. Three of the four features of the present MHD 
schemes summarized in the abstract are validated by typical 1-D and 2-D MHD test cases. 
The new method of defining high-resolution shock-capturing schemes for the conservative 
MHD equations exhibits smaller V-B numerical error when compared to the standard 
numerical method for solving the non-conservative MHD equations. The approach of using 
the non- conservative eigensystem when solving the conservative equations is also applicable 
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WAV66+AD8 101x101 



WAV66+AD8 401x 401 



WAV66+AD8 201x 201 



WAV66+AD8 801x 801 



Fig. 21. Temperature contours of the vortex pairing problem. Grid refinement for the 
WAV66+AD8 scheme. Conservative MHD equations. 
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Fig. 22. Vortex pairing problem. V B vs. time. Harten-Yee (solid) WAV66+AD8 (dashed). 

in the context of commonly used shock-capturing schemes for the MHD equations in the 
literature. 

More importantly, the filter approach of our new method provides a natural and effi- 
cient way to minimize the V- B numerical error. It also offers a variety of options to cater 
to different flow types of the problem. Since there are many variants in the minimization 
of the divergence of the magnetic field (V*B) numerical error, Part II of the companion 
paper [21] is devoted to extensive comparison of these variants. In yet another forthcoming 
paper, Part III, complex multiscale high speed turbulent astrophysical applications will 
be sought. 
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